#!/usr/bin/perl -w
#use strict;
#use warnings;

## VERSION 3.01 ## ajout d'une fonction de tri ## JCC avril 2013 ##

my ($USAGE) =  "\n\n\t****** USAGE of preprocess.pl PROGRAM ******\n\n\n\n\tUSAGE: perl $0 <panam.ini file> \n\n\n\n";
die "$USAGE" if ((scalar(@ARGV))< 1);
my $option_file = $ARGV[0];
chomp($option_file);

die "\n\n\t=> Cannot find configuration file: $option_file.\n\n" unless (-e $option_file);
die "\n\n\t=> Configuration file, $option_file, appears to be empty!\n\n" if (-z $option_file);
open(OFILE, "<$option_file") || die "Cannot open input file : $option_file\n";
# 
my $usearch;
my $seqIn;
my $barIn;
my $Clst;
my $NGS_id_Results;
my $preprocess_output = "preprocess_output";
my $barcode;
my $pool;
my $us_version = "3.0.617";
my $path = "/home/panam/panam_v3/bin/uclust3.0.617";

my $min_norm; 
my $nbr_seq_norm; 

my $lengthCut;
my $maxlength ;
########################"""

my $nbrBar; my $ficBar;

while (<OFILE>) {
	chomp;
	my $line = $_;

# discard blank line
	if ($line =~ m/^\s*$/g) {
		next;
	}

# discard comment line
	elsif($line =~ m/^\s*#/g) {
		next;
	}

# get input options line
        else {
	
		if($line=~m/SEQUENCE_FILE/gi) {
			$line=~s/SEQUENCE_FILE//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) {
				$seqIn = $line; #print "$inputSeqName\n";
			}
			chomp ($seqIn);
			# See if the 454 Sequence file exist and contains something
			die "\n\n\tCannot find 454 sequence file: $seqIn. Check panam.ini.\n\n" unless (-e $seqIn);
			die "\n\n\t454 sequence file, $seqIn, appears to be empty!\n\n" if (-z $seqIn);

			my $ligne1= `sed -n '1p' $seqIn`;
			my $ligne2= `sed -n '2p' $seqIn`;	
# 			if (($ligne1 !~ /^>/) or ($ligne2 !~ /^[ATCGU]/)) { #### modif 25/5
			if (($ligne1 !~ /^>/) or ($ligne2 !~ /^[ATCGU]/i)) {
				die "\n\n\tYour file $seqIn does not seem to contain sequences in fasta format.\n\n";
			}
		}

		if($line=~m/BARCODE_FILE/gi) { 
			$line=~s/BARCODE_FILE//gi; ;
			$line=~s/\t//gi; 
			$line=~s/^\s+//gi; 
			$line=~s/\s+$//gi; 
			if($line) {
				$barIn = $line; #print "$barIn\n";
			}

	
			if (!(defined $barIn ) or ($barIn eq "")) {
				$barcode =0;
				$nbrBar = 0;
				########ajout 13/5
				$ficBar = "no barcode file" ;
				#############fin ajout 13/5
			}
			elsif ((defined $barIn) and ($barIn ne "")){
				$barcode = 1;
				##########" ajout 13/5
				$ficBar = $barIn;
				########fin ajout 13/5     

				# See if the barcode file exist and contains something
				die "\n\n\tCannot find barcode file: $barIn. Check panam.ini.\n\n" unless (-e $barIn);
				die "\n\n\tBarcode file, $barIn, appears to be empty!\n\n" if (-z $barIn);

				open (BB, $barIn);
				my @barIN= <BB>;
				close BB;
				foreach my $e (@barIN) {
					if ($e !~ /\t/) {die "\n\n\tYour bar code file $barIn does not seem to be in the right format.\n\n";}
					my @ty = split(/\t/, $e);
					if (($ty[0] =~ /^\d/) or ($ty[0]=~ /_/) or ($ty[0]=~ /\//)) { die "\n\n\tThe bar code Ids in $barIn should start with no number and contain no underscorses.\n\n";}
					if ($ty[1] !~ /[ATCGU]/i) {die "\n\n\tYour bar code file $barIn does not seem to be in the right format.\n\n";}
				}
				############### fin modif + ajout

				my $nbr = `wc -l $barIn`;
				if($nbr =~ /(\d.?)\s/) { $nbrBar = $1} #### modif 25/5 w --> d 
			}
		}

		################
		if($line=~m/LOWER_SEQUENCE_LENGTH_CUTOF/gi) {
			$line=~s/LOWER_SEQUENCE_LENGTH_CUTOF//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) {
				$lengthCut = $line; #print "$lengthCut\n";
			}
			die "\n\n\t Length cutoff missed. Check panam.ini.\n\n" unless ($lengthCut ne "");
		}
 
		if($line=~m/MAX_SEQUENCE_LENGTH/gi) {
			$line=~s/MAX_SEQUENCE_LENGTH//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) {
				$maxlength = $line; #print "$maxlength\n";
			}
			die "\n\n\t Length cutoff missed. Check panam.ini.\n\n" unless ($maxlength ne "");
		}

		############### 
	
		if($line=~m/CLUSTERING_CUTOFF/gi) {
			$line=~s/CLUSTERING_CUTOFF//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) {
				$Clst = $line; #print "$Clst\n";
			}
		die "\n\n\t Clustering cutoff missed. Check panam.ini.\n\n" unless ($Clst <=1);
		die "\n\n\t Invalid input! clustering cutoff must be specified as a fractional identity in the range 0.0 to 1.0\n\n" unless (($Clst ne "") and ($Clst >=0));
		}
	
		if($line=~m/454_RUN_IDENTIFIER/gi) { ###### dossier !
			$line=~s/454_RUN_IDENTIFIER//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) {
				$NGS_id_Results = $line; #print "$NGS_id_Results\n";
			}
	 		die "\nOutput directory for NGS analyses is not defined. Check panam.ini.\n\n" unless ( defined $NGS_id_Results);
			unless (-d $NGS_id_Results) { mkdir $NGS_id_Results || die "Could not create Directory $NGS_id_Results $!\n"; }
		}

		if (-d $NGS_id_Results."/".$preprocess_output) { `rm -r $NGS_id_Results/$preprocess_output` } 
		if (!(-d $NGS_id_Results."/".$preprocess_output)) { mkdir $NGS_id_Results."/".$preprocess_output || die "\n\tCould not create Directory $NGS_id_Results/$preprocess_output $!\n"; }

		if($line=~m/POOLING/gi) {
			$line=~s/POOLING//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) {
				$pool = $line; 
			}
			die "\n\n\t pooling value missed. Check panam.ini.\n\n" unless (($pool ne ""));
			die "\n\n\t Invalid pooling value!\n\n" unless (($pool eq "yes") or ($pool eq "no"));	
		}

		unless (-d $NGS_id_Results."/".$preprocess_output) { mkdir $NGS_id_Results."/".$preprocess_output || die "Could not create Directory $NGS_id_Results/$preprocess_output $!\n"; }

		############### normalisation

		if($line=~m/NBR_SEQ_NORM/gi) {
			$line=~s/NBR_SEQ_NORM//gi;
			$line=~s/\t//gi;
			$line=~s/^\s+//gi;
			$line=~s/\s+$//gi;
			if($line) {
				$nbr_seq_norm = $line;
			} 
		}
	}
}

die "\n\n\t Invalid number of sequences to pick for the normalization! Check panam.ini.\n\n" if ((defined $nbr_seq_norm) and ($nbr_seq_norm !~ /^\d+$/));
if (!(defined $nbr_seq_norm)) { $nbr_seq_norm = '200'}

##################### Ajout JCC avril 2013 : by sample order

my $socat = `cat -n $barIn | awk -F"\t" '{ print\$1"\t"\$2 }'`;
my @socat = split("\n", $socat);
my %so = ();

foreach my $line (@socat) {
    my ($val, $key) = split("\t", $line);
    $so{$key} = $val;
}


###################  barcoding
my $dirName = $NGS_id_Results."/".$preprocess_output ;

if ($barcode == 0 ) {
	open (SI, $seqIn) || die "can not open file";
	my %seq; my $a;

	while (<SI>) {
	   if ($_=~ />(.*?)\s+/) {  
	      $a = $1
	   }
	    else {
		$seq{$a}.=$_
	   }
	}

	open (SO, ">".$dirName."/sample_sample.fasta"); 
	foreach my $k (keys %seq) {
	  print SO ">sample_$k\n$seq{$k}\n";
	}
	close SO;

	#$barIn{'sample'}=1;
}

my %sample_name;
if ($barcode == 1) {
	my $l;
	my %barcode;

	open (B, $barIn) || die "can not open file"; 

	while (<B>) {
		my $line = $_;
		$line=~s/\s+$//gi;
		my @t = split (/\t/, $line); 
		chomp ($t[0]); chomp ($t[1]); 
		$barcode{$t[0]} = $t[1];
		$l = length($t[1]); 
	}
	close B;

	my %seq; my $aa;
	open(S, $seqIn) || die "can not open file $seqIn";
	while (<S>) {
		if ($_=~ />(.*?)\n/) {
			$aa=$1 ; 
		}
		else {
			chomp($_);
			$seq{$aa}.=$_
		}
	}
	close S;

	my $i=0;
	my %posF; my %posR; 
	
	my %pos;
	my $startF; my $startR;	

	open (FF, $NGS_id_Results."/tmp/outputF_Ids_tmp.fuzznuc") || die "can not open file" ;
	while (<FF>) {
		if ($_=~ /SeqName/) { }
		else {
			my @t= split (/\t/,$_); 
			$pos{$t[0]}{'F'}{'start'} = $t[2]; 
			$pos{$t[0]}{'F'}{'end'} = $t[3];
 
			##### récupérer les positions de startF et startR pour les comparer, connaître le primer qui démarre le séquençage et chercher le tag en amont de cette amorce. 
			if (exists $pos{$t[0]}{'F'}{'start'}) {$startF = $pos{$t[0]}{'F'}{'start'}}	
		}
	}
	close FF;
	

	open (FR, $NGS_id_Results."/tmp/outputR_Ids_tmp.fuzznuc") || die "can not open file" ;
	while (<FR>) {
		if ($_=~ /SeqName/) { }
		else {
			my @t= split (/\t/,$_);
			$pos{$t[0]}{'R'}{'start'} = $t[2];
			$pos{$t[0]}{'R'}{'end'} = $t[3];
			if (exists $pos{$t[0]}{'R'}{'start'}) {$startR = $pos{$t[0]}{'R'}{'start'}}	
		}
	}
	close FR;


#############################################

	my %amont; my %res; my %posAbs;
	foreach my $k (keys %seq) { 
		my $m; my $f;
		if ((defined $pos{$k}{'F'}{'start'}) and (defined $pos{$k}{'R'}{'end'})) {
			if ($startF < $startR) {	
				$m = $pos{$k}{'F'}{'start'} -1;
				$f = $pos{$k}{'R'}{'end'}; 
			}
			else {
				$m = $pos{$k}{'R'}{'start'} -1;
				$f = $pos{$k}{'F'}{'end'}
			}
		}
		elsif ((defined $pos{$k}{'F'}{'start'}) and (!(defined $pos{$k}{'R'}{'end'}))) {
			$m = $pos{$k}{'F'}{'start'} -1;
		}
		#else {
		elsif ((!(defined $pos{$k}{'F'}{'start'})) and (defined $pos{$k}{'R'}{'end'})) {
			$m = $pos{$k}{'R'}{'start'} -1;
		}
	
		$posAbs{$k}{'start'} = $m;
		if (defined $f) {	
			$posAbs{$k}{'end'} = $f;	
			#print "m --- $m\nf --- $f\n";
		}

		$amont{$k}= substr($seq{$k}, $m - $l, $l); 
		foreach my $bar (keys %barcode) { 
			if ($barcode{$bar} eq $amont{$k}) {  
				$res{$bar}{$k} = $seq{$k}; 
			}
		}
	}	

	my %clean_seq;

	#######################################

	foreach my $k (keys %res) {
		foreach my $p (keys %{$res{$k}}) {
			my $clean_seq;
			my $length;
			if (exists $posAbs{$p}{'end'}) {
				$length = $posAbs{$p}{'end'} - $posAbs{$p}{'start'};
			}
			else {
				$length = $maxlength
			}
		
			$clean_seq = substr ($seq{$p}, $posAbs{$p}{'start'}, $length);
			my $ll = length ($clean_seq);
			if ($ll >= $lengthCut) {
				$clean_seq{$k}{$p}=$clean_seq;
			}
		}
	}
	
	#########################################"
	
	foreach my $e (keys %clean_seq) { 
		open ($e, ">$dirName/"."sample_".$e.".fasta") or die $!;
		foreach my $k (keys %{$clean_seq{$e}}) {
			my $l = length ($clean_seq{$e}{$k});
			print $e ">".$e."_"."$k\n$clean_seq{$e}{$k}\n"
		}
		close $e;
	}
}


$DIR = $dirName;

#######"""check if $DIR is empty 
my $dir = `ls -A $DIR | wc -c ` ;
die "\n$DIR seems to be empty.\n\n" unless ( $dir != 0);

# ############################### Clustering and diversity index

opendir DIR, $DIR || die "can not open $DIR";

my @files = readdir DIR;
my $sca = scalar (@files); 
my %index; my %nb_clust;

print "Clustering and rarefaction ...\n";
my %seq;
$min_norm=$nbr_seq_norm ;
foreach my $f (@files) {
	if (($f=~ /sample_(.*?).fasta/) and ($f ne "nobarcode.fasta")) {
		my %seq_sample;
		my $e = $1;
		print "Sample ---- $e\n";

		########### ajout 3/4/2012
		#$sample_name{$e}=1; ### modif 26/7
		########## fin ajout

		if ( !(-z $DIR."/sample_".$e.".fasta")) {  ### modif 26/7
		$sample_name{$e}=1;  ### modif 26/7

		open (S, $DIR."/sample_".$e.".fasta") || die "can not open file";
		my $aa;
		while (<S>) {
			if ($_=~ /^>(.*?)\n/) {
				$aa = $1
			}

			else {
				chomp($_);
				$seq{$aa}.=$_;
				########### ajout 29/3/2012
				$seq_sample{$aa}.=$_
				########## fin ajout
			}
		}
		close S;

		################ 
		my $keys_seq_sample = keys (%seq_sample);
		$index{$e}{'nbrseq'} = $keys_seq_sample; #print "$e --- $index{$e}{'nbrseq'}\n";
		################ 

		if (($us_version eq "1.1.579q") or ($us_version eq "3.0.617" )) {
			`$path --sort $DIR"/sample_"$e".fasta" --output $DIR"/"$e"_sorted" 2>&1`;
			`$path --input $DIR"/"$e"_sorted" --uc $DIR"/"$e.uc --id $Clst --optimal --rev 2>&1`;
		}
		
		else  {
			`$path --sort $DIR"/sample_"$e".fasta" --output $DIR"/"$e"_sorted" 2>&1` ;
			`$path --cluster $DIR"/"$e"_sorted" --id $Clst --uc $DIR"/"$e.uc --optimal --rev 2>&1`
		}

		my $resultUC = $DIR."/".$e.".uc";		

		die "\n\n\tCannot find USEARCH output: $resultUC.\n\n" unless (-e $resultUC);
		die "\n\n\tUSEARCH output, $resultUC, appears to be empty!\n\n" if (-z $resultUC);
		open (F, $DIR."/".$e.".uc") ;

		my %seed; my %seed_seq;
		my $n=0; my %read_to_cluster; my $max=0;

		while (<F>) {
			my @t = split (/\t/, $_);

			if ($t[0] eq "H") {
				$n++;
				$read_to_cluster{$e}{$n}=$t[1];
				chomp ($t[9]); chomp ($t[8]);
				$seed_seq{$t[9]} .= $t[8].", " ;
				$seed{$t[9]}++ ;
			}

			if ($t[0] eq "S") {
				$n++;
				$read_to_cluster{$e}{$n}=$t[1];	
				$seed_seq{$t[8]} .= $t[8].", " ;	
				$seed{$t[8]}++ ;
			}
		}

		$index{$e}{'nbrOTU'} = keys (%seed);

		################## generate seeds 'files (*_Seq_OTU) and nbr_seq / OTU (*_OTU)

		open (R, ">".$DIR."/".$e."_OTU");
		open (SEQ, ">".$DIR."/".$e."_Seq_OTU");
		my $i=1; my $n1=0; my $n2=0; 	
		my $ninfeg10=0; my $nsup10=0; my $Nrare=0; my $gamma;
		my $n3=0; my $n4=0; my $n5=0; my $n6=0; my $n7=0; my $n8=0; my $n9=0; my $n10=0;
		print R "OTU\tSeed sequence\t#Sequences\n";

		foreach my $e ((sort { $seed {$b} <=> $seed {$a}} keys %seed) ){
			if ($seed{$e} == 1 ) {$n1++}
			if ($seed{$e} == 2 ) { $n2++}

		############### ace

			if ($seed{$e} == 3 ) {$n3++}
			if ($seed{$e} == 4 ) {$n4++}
			if ($seed{$e} == 5 ) {$n5++}	
			if ($seed{$e} == 6 ) {$n6++}	
			if ($seed{$e} == 7 ) {$n7++}
			if ($seed{$e} == 8 ) {$n8++}	
			if ($seed{$e} == 9 ) {$n9++}
			if ($seed{$e} == 10 ) {$n10++}

			if ($seed{$e} <= 10) {$ninfeg10++}	
			if ($seed{$e} > 10) {$nsup10++}		

			print R "$i\t$e\t$seed{$e}\t$seed_seq{$e}\n";
			$i++;
			print SEQ ">$e\n$seq{$e}\n";
# 			print OTU ">$e\n$seq{$e}\n";
		}

		close R; 
		close SEQ;

		######################rarefaction

		my $seuil = $index{$e}{'nbrseq'}/10;
		#my $int = int ($seuil);

		#################
		my $int;
		if ($seuil >= 1) {
			$int = int ($seuil)
		}
		else {
			$int = 1
		} 
		#######################

		if ($n> $max) {
			$max = $n
		}

		my $nb_a_ech=0;
		my $nb_ech=int($n/$int); 

		while($nb_ech>0){
			$nb_a_ech+=$int;

			my %vu_clust;
			my %vu_read;
			my $reste_a_ech=$nb_a_ech;

			while($reste_a_ech>0){
				my $ech=rand()*$n;
				$ech=int($ech);
				if (defined($vu_read{$e}{$ech})){}
				else{
					$vu_read{$e}{$ech}=1;
					if( defined $read_to_cluster{$e}{$ech}) {
						 $vu_clust{$read_to_cluster{$e}{$ech}}++;
					}
					$reste_a_ech--;
				}
			}
			my @liste_clust=keys %vu_clust;
			$nb_clust{$e}{$nb_a_ech}=$#liste_clust+1;
			$nb_ech--;
		}

		##################Schao

		if ((($n1 >0) and ($n2 >= 0)) or (($n1==0 ) and ($n2 == 0))) {
			$index{$e}{'chao'} = $index{$e}{'nbrOTU'} + (($n1*($n1 - 1))/(2*($n2+1)));
		}

		else {
			$index{$e}{'chao'} = $index{$e}{'nbrOTU'} + (($n1*$n1)/(2*$n2));
		}

		$index{$e}{'chao'} = sprintf("%.2f", $index{$e}{'chao'});

		###############coverage

	#	my $cov = 1- ($n1/$index{$e}{'nbrseq'});
	#	$cov = $cov * 100;
	#	$index{$e}{'cov'} = sprintf("%.2f", $cov);

		my $cov;
		if ($index{$e}{'nbrseq'} != 0) {
			$cov = 1- ($n1/$index{$e}{'nbrseq'});
			$cov = $cov * 100;
			$index{$e}{'cov'} = sprintf("%.2f", $cov);
		}
		else {
			$index{$e}{'cov'} = "NA"
		}

		##################Sace

		my @rare = ($n1, $n2, $n3, $n4, $n5, $n6, $n7, $n8, $n9, $n10);

		for (my $i=0; $i< 10; $i++ ) {
			my $u = ($i+1)*$rare[$i];
			$Nrare += $u
		}
		for (my $i=0; $i< 10; $i++ ) {
			my $v = ($i+1)*($i+1-1)*$rare[$i] ;
			$gamma += $v
		}
			
		#my $Cace = 1- $n1/$Nrare;

		#my $bst = $ninfeg10*$gamma; my $mkm = $Cace*$Nrare*($Nrare-1);
		#my $Gace = ($bst/$mkm) - 1;
		
		#if ($Gace < 0) {$Gace = 0}		

		#my $ace  = $nsup10 + ($ninfeg10/$Cace) + ($n1/$Cace)*$Gace ;
		#$index{$e}{'Sace'} = sprintf("%.2f", $ace);

		################################# modif
		
		my $Cace; my $mkm;
		if ($Nrare != 0) {
			$Cace = 1 - $n1/$Nrare;
			my $bst = $ninfeg10*$gamma; 
			$mkm = $Cace*$Nrare*($Nrare-1);
			my $Gace; 
			if ( $mkm != 0) {
				$Gace = ($bst/$mkm) - 1;
				if ($Gace < 0) {$Gace = 0}
			}
			my $ace; 
			if ($Cace != 0) {
				$ace = $nsup10 + ($ninfeg10/$Cace) + ($n1/$Cace)*$Gace ;
				$index{$e}{'Sace'} = sprintf("%.2f", $ace);
			}
		}
		if (($Nrare ==0) or($mkm ==0) or ($Cace ==0)) {
			$index{$e}{'Sace'} = "NA";
		}	
	
		
# 		##########################shannon

		open (SH, $DIR."/".$e."_OTU") ;
		my $sum = 0; my $sum2 = 0; my $sh;
		my @tab = <SH>;

		foreach my $elt (@tab) {
			if ($elt !~ /OTU/) {
				my @t = split (/\t/, $elt) ;
				chomp ($t[2]);
				$sum += $t[2];
			}
		}

		#foreach my $el (@tab) {
		#	if ($el !~ /OTU/) { 
		#		my @t = split (/\t/, $el);
		#		chomp ($t[2]);
		#		$t[2] = $t[2]/$sum;
		#		$t[2] = $t[2] * log($t[2]);
		#		$sh += $t[2];
		#	}
		#}

		#$sh = $sh*(-1);
		#$index{$e}{'shannon'} = sprintf("%.2f", $sh);
		#close SH;

		if ($sum != 0) {
			foreach my $el (@tab) {
				if ($el !~ /OTU/) { 
					my @t = split (/\t/, $el);
					chomp ($t[2]);
					$t[2] = $t[2]/$sum;
					$t[2] = $t[2] * log($t[2]);
					$sh += $t[2];
				}
			}

			$sh = $sh*(-1);
			$index{$e}{'shannon'} = sprintf("%.2f", $sh);
		}
		else {
			$index{$e}{'shannon'} = "NA"
		}

		close SH;
	}
	}
}

print "Clustering finished.\nPrinting results ...\n";

if (($barcode == 0) or ($nbrBar <=1)) {
	open (IN, ">".$NGS_id_Results."/index_".$Clst) ;
	print IN "Sample\t#Sequences\t#OTUs\tChao1\tShannon_Index\tACE\tCoverage\n"; ##### modif 31/5
  
#	foreach my $e (sort keys %index) {
	foreach my $e (sort {$so{$a} <=> $so{$b}} keys %index) { ## ajout d'une fonction de tri JCC avril 2013
		open (RAR, ">".$DIR."/".$e."_rarefaction");
		print RAR "#Sequences\t#OTUs\n";
		print IN "sample\t$index{$e}{'nbrseq'}\t$index{$e}{'nbrOTU'}\t$index{$e}{chao}\t$index{$e}{'shannon'}\t$index{$e}{'Sace'}\t$index{$e}{'cov'}\n";
		foreach my $p (sort {$a <=> $b} keys %{$nb_clust{$e}})  {
			print RAR "$p\t$nb_clust{$e}{$p}\n";
		}
		close RAR;
	} 
	close IN;
}

else {
############################### normalisation
	print "Normalization ...\n";
	my %index_norm; my %nb_clust_norm ;
	
	my $DIR_norm = $DIR."/norm";
	`mkdir $DIR_norm`;

	foreach my $e (keys %sample_name) {
		my %seq_sample;

		if (!(-z $DIR."/sample_".$e.".fasta")) {
		open (S, $DIR."/sample_".$e.".fasta") || die "can not open file";
		my $aa;
		while (<S>) {
			if ($_=~ /^>(.*?)\n/) {
				$aa = $1
			}
			else {
				chomp($_);
				$seq_sample{$aa}.=$_;
			}
		}
		close S;

	
		my $keys_seq_sample = keys (%seq_sample);  

		if ($keys_seq_sample >= $min_norm) {       
			my $norm=0;
			open (SS, ">".$DIR_norm."/sample_".$e."_".$min_norm.".fasta") || die "can not open file";
			foreach my $k (keys %seq_sample) {
				print SS ">$k\n$seq_sample{$k}\n";
				$norm ++;
				if ($norm == $min_norm) {
					last
				}
			}
			close SS;
			
			$index_norm{$e}{'nbrseq'} = $min_norm;

			if (($us_version eq "1.1.579q") or ($us_version eq "3.0.617" )) {
				`$path --sort $DIR_norm"/sample_"$e"_"$min_norm".fasta" --output $DIR_norm"/"$e"_"$min_norm"_sorted" 2>&1`;
				`$path --input $DIR_norm"/"$e"_"$min_norm"_sorted" --uc $DIR_norm"/"$e"_"$min_norm".uc" --id $Clst --optimal --rev 2>&1`;
			}
				
			else  {
				`$path --sort $DIR_norm"/sample_"$e"_"$min_norm".fasta" --output $DIR_norm"/"$e"_"$min_norm"_sorted" 2>&1` ;
				`$path --cluster $DIR_norm"/"$e"_"$min_norm"_sorted" --id $Clst --uc $DIR_norm"/"$e"_"$min_norm".uc" --optimal --rev 2>&1`
			}

			my $resultUC = $DIR_norm."/".$e."_".$min_norm.".uc";		
			
			die "\n\n\tCannot find USEARCH output: $resultUC.\n\n" unless (-e $resultUC);
			die "\n\n\tUSEARCH output, $resultUC, appears to be empty!\n\n" if (-z $resultUC);
			open (FNORM, $DIR_norm."/".$e."_".$min_norm.".uc") ;

			my %seed; my %seed_seq;
			my $n=0; my %read_to_cluster; my $max=0;

			while (<FNORM>) {
				my @t = split (/\t/, $_);
					if ($t[0] eq "H") {
					$n++;
					$read_to_cluster{$e}{$n}=$t[1];
					chomp ($t[9]); chomp ($t[8]);
					$seed_seq{$t[9]} .= $t[8].", " ;
					$seed{$t[9]}++ ;
				}

				if ($t[0] eq "S") {
					$n++;
					$read_to_cluster{$e}{$n}=$t[1];	
					$seed_seq{$t[8]} .= $t[8].", " ;	
					$seed{$t[8]}++ ;
				}
			}

			$index_norm{$e}{'nbrOTU'} = keys (%seed);
			
			################## generate seeds 'files (*_Seq_OTU) and nbr_seq / OTU (*_OTU)

			open (RNORM, ">".$DIR_norm."/".$e."_OTU"); #### modif 4/6
			open (SEQNORM, ">".$DIR_norm."/".$e."_"."$min_norm"."_Seq_OTU");

			my $i=1; my $n1=0; my $n2=0; 	
			my $ninfeg10=0; my $nsup10=0; my $Nrare=0; my $gamma;
			my $n3=0; my $n4=0; my $n5=0; my $n6=0; my $n7=0; my $n8=0; my $n9=0; my $n10=0;

			print RNORM "OTU\tSeed sequence\t#Sequences\n";

			foreach my $e ((sort { $seed {$b} <=> $seed {$a}} keys %seed) ){
				if ($seed{$e} == 1 ) {$n1++}
				if ($seed{$e} == 2 ) { $n2++}

			############### ace

				if ($seed{$e} == 3 ) {$n3++}
				if ($seed{$e} == 4 ) {$n4++}
				if ($seed{$e} == 5 ) {$n5++}	
				if ($seed{$e} == 6 ) {$n6++}	
				if ($seed{$e} == 7 ) {$n7++}
				if ($seed{$e} == 8 ) {$n8++}	
				if ($seed{$e} == 9 ) {$n9++}
				if ($seed{$e} == 10 ) {$n10++}

				if ($seed{$e} <= 10) {$ninfeg10++}	
				if ($seed{$e} > 10) {$nsup10++}		

					print RNORM "$i\t$e\t$seed{$e}\t$seed_seq{$e}\n";
					$i++;
					print SEQNORM ">$e\n$seq{$e}\n";

				}

			close RNORM; 
			close SEQNORM;
# 
			######################rarefaction

			my $seuil = $index_norm{$e}{'nbrseq'}/10;
			#my $int = int ($seuil);

			#################
			my $int;
			if ($seuil >= 1) {
				$int = int ($seuil)
			}
			else {
				$int = 1
			} 
			#######################

			if ($n> $max) {
				$max = $n
			}

			my $nb_a_ech=0;
			my $nb_ech=int($n/$int); 

			while($nb_ech>0){
				$nb_a_ech+=$int;
				my %vu_clust;
				my %vu_read;
				my $reste_a_ech=$nb_a_ech;

				while($reste_a_ech>0){
					my $ech=rand()*$n;
					$ech=int($ech);
					if (defined($vu_read{$e}{$ech})){}
					else{
						$vu_read{$e}{$ech}=1;
						if( defined $read_to_cluster{$e}{$ech}) {
							$vu_clust{$read_to_cluster{$e}{$ech}}++;
						}
						$reste_a_ech--;
					}
				}
				my @liste_clust=keys %vu_clust;
				$nb_clust_norm{$e}{$nb_a_ech}=$#liste_clust+1;
				$nb_ech--;
			}

			##################Schao

			if ((($n1 >0) and ($n2 >= 0)) or (($n1==0 ) and ($n2 == 0))) {
				$index_norm{$e}{'chao'} = $index_norm{$e}{'nbrOTU'} + (($n1*($n1 - 1))/(2*($n2+1)));
			}

			else {
				$index_norm{$e}{'chao'} = $index_norm{$e}{'nbrOTU'} + (($n1*$n1)/(2*$n2));
			}

			$index_norm{$e}{'chao'} = sprintf("%.2f", $index_norm{$e}{'chao'});

			###############coverage

		#	my $cov = 1- ($n1/$index_norm{$e}{'nbrseq'});

			my $cov;
	  		if ($index_norm{$e}{'nbrseq'} != 0) {
			  $cov = 1- ($n1/$index_norm{$e}{'nbrseq'});
			  $cov = $cov * 100;
			  $index_norm{$e}{'cov'} = sprintf("%.2f", $cov);
			 }
	 		else {
				$index_norm{$e}{'cov'} = "NA"
	 		}	

			##################Sace

			my @rare = ($n1, $n2, $n3, $n4, $n5, $n6, $n7, $n8, $n9, $n10);
			for (my $i=0; $i< 10; $i++ ) {
				my $u = ($i+1)*$rare[$i];
				$Nrare += $u
			}
			for (my $i=0; $i< 10; $i++ ) {
				my $v = ($i+1)*($i+1-1)*$rare[$i] ;
				$gamma += $v
			}
			
			#my $Cace = 1- $n1/$Nrare;

			#my $bst = $ninfeg10*$gamma; my $mkm = $Cace*$Nrare*($Nrare-1);
			#my $Gace = ($bst/$mkm) - 1;
		#	
			#if ($Gace < 0) {$Gace = 0}		
			#my $ace  = $nsup10 + ($ninfeg10/$Cace) + ($n1/$Cace)*$Gace ;
			#$index_norm{$e}{'Sace'} = sprintf("%.2f", $ace);
	
			my $Cace; my $mkm; 
			if ($Nrare != 0) {
				$Cace = 1 - $n1/$Nrare;
				my $bst = $ninfeg10*$gamma; 
				$mkm = $Cace*$Nrare*($Nrare-1);
				my $Gace; 
				if ( $mkm != 0) {
					$Gace = ($bst/$mkm) - 1;
					if ($Gace < 0) {$Gace = 0}
				}
				my $ace; 
				if ($Cace != 0) {
					$ace = $nsup10 + ($ninfeg10/$Cace) + ($n1/$Cace)*$Gace ;
					$index_norm{$e}{'Sace'} = sprintf("%.2f", $ace);
				}
			}
			if (($Nrare ==0) or ($mkm ==0) or ($Cace == 0)) {
				$index_norm{$e}{'Sace'} = "NA";
			}
	
		# 	##########################shannon

			open (SHNORM, $DIR_norm."/".$e."_OTU") ;
			my $sum = 0; my $sum2 = 0; my $sh;
			my @tab = <SHNORM>;

			foreach my $elt (@tab) {
				if ($elt !~ /OTU/) {
					my @t = split (/\t/, $elt) ;
					chomp ($t[2]);
					$sum += $t[2];
				}
			}

			if ($sum != 0) {
				foreach my $el (@tab) {
			 		 if ($el !~ /OTU/) { 
						  my @t = split (/\t/, $el);
						  chomp ($t[2]);
						  $t[2] = $t[2]/$sum;
						  $t[2] = $t[2] * log($t[2]);
						  $sh += $t[2];
					  }
		 		 }

		 		 $sh = $sh*(-1);
		 		 $index_norm{$e}{'shannon'} = sprintf("%.2f", $sh);
	  		}

	 		else {
				$index_norm{$e}{'shannon'} = "NA"
			}

			close SHNORM;
		}
		}
	}     

	open (IN, ">".$NGS_id_Results."/index_".$Clst) ;
	print IN "Sample\t#Sequences\tNormalized_Sequences\t#OTUs\tNormalized_OTU\tChao1\tNormalized_Chao1\tShannon_Index\t Normalized_Shannon_Index\tACE\tNormalized_ACE\tCoverage\tNormalized_Coverage\n"; 

	foreach my $e (sort {$so{$a} <=> $so{$b}} keys %index) { ## VERSION 3.01 ##
		open (RAR, ">".$DIR."/".$e."_rarefaction");
		print RAR "#Sequences\t#OTUs\n";

		open (RARNORM, ">".$DIR_norm."/".$e."_".$min_norm."_rarefaction");
		print RARNORM "#Sequences\t#OTUs\n";

		if (defined $index_norm{$e}) {
			print IN "$e\t$index{$e}{'nbrseq'}\t$index_norm{$e}{'nbrseq'}\t$index{$e}{'nbrOTU'}\t$index_norm{$e}{'nbrOTU'}\t$index{$e}{chao}\t$index_norm{$e}{chao}\t$index{$e}{'shannon'}\t$index_norm{$e}{'shannon'}\t$index{$e}{'Sace'}\t$index_norm{$e}{'Sace'}\t$index{$e}{'cov'}\t$index_norm{$e}{'cov'}\n";
		}
		else {
			print IN "$e\t$index{$e}{'nbrseq'}\tNA\t$index{$e}{'nbrOTU'}\tNA\t$index{$e}{chao}\tNA\t$index{$e}{'shannon'}\tNA\t$index{$e}{'Sace'}\tNA\t$index{$e}{'cov'}\tNA\n";
		}

		foreach my $p (sort {$a <=> $b} keys %{$nb_clust{$e}})  {
			print RAR "$p\t$nb_clust{$e}{$p}\n";	
		}
		close RAR;

		
		foreach my $pp (sort {$a <=> $b} keys %{$nb_clust_norm{$e}})  {
			print RARNORM "$pp\t$nb_clust_norm{$e}{$pp}\n";	
		}
		close RARNORM;
	}

	close IN;
}

#########################""""""fin ajout 3/4/2012

######################################"
open (L ,">>$NGS_id_Results/PREPROCESS_LOG");

######"remplacé $barIN par $ficBar 13/5
print L "Barcode file: $ficBar\nClustering threshold: $Clst\n";

if ($pool eq "no") {
 	my $dir = $NGS_id_Results."/".$preprocess_output ; ##### modif 3/6

	opendir DIR, $dir || die "can not open $dir";
	my @files = readdir DIR ;
	foreach my $f (@files) { 
		if ($f =~ /Seq_OTU/) {
# 			`cat $NGS_id_Results/$preprocess_output/$f >>  $NGS_id_Results/All_Seq_OTU_$Clst` ; ##### modif 3/6
			 `cat $dir/$f >>  $NGS_id_Results/All_Seq_OTU_$Clst` 
		}
	}

	print "\nSamples description printed in index file.\nOTUs printed in All_Seq_OTU_$Clst file.\nYou can now process All_Seq_OTU_$Clst with panam.pl\n\n";
}

########################## pooled_sample ##########################

if ($pool eq "yes") {
	print "Processing pooled samples ...\n";
	`mkdir $NGS_id_Results/$preprocess_output/pooled_sample`;
	my $pooled_sample = $NGS_id_Results."/".$preprocess_output."/pooled_sample";
	
	foreach my $f (@files) {
		if (($f ne "nobarcode.fasta") and ($f =~ /sample_.*?fasta/)) {
 			`cat $NGS_id_Results/$preprocess_output/$f >> $pooled_sample/pooled_sample_Seq.fasta` ; 
#			`cat $dir/$f >> $pooled_sample/pooled_sample_Seq.fasta` ;
		}
	}


	if (($us_version eq "1.1.579q") or ($us_version eq "3.0.617" )) {
		`$path --sort $pooled_sample/pooled_sample_Seq.fasta --output $pooled_sample/pooled_sample_Seq_sorted.fasta 2>&1`;
		`$path --input $pooled_sample/pooled_sample_Seq_sorted.fasta --uc $pooled_sample"/pooled_sample_Seq_"$Clst".uc" --id $Clst --optimal --rev 2>&1`;
	}
		
	else  {
		`$path --sort $pooled_sample/pooled_sample_Seq.fasta --output $pooled_sample/pooled_sample_Seq_sorted.fasta 2>&1` ;
		`$path --cluster $pooled_sample/pooled_sample_Seq_sorted.fasta --id $Clst --uc $pooled_sample"/pooled_sample_Seq_"$Clst".uc" --optimal --rev 2>&1`
	}
	
	open (RPS, ">".$pooled_sample."/pooled_sample_OTU");
	open (SEQPS, ">".$pooled_sample."/pooled_sample_SEQ_OTU");
	open (FPS, $pooled_sample."/pooled_sample_Seq_".$Clst.".uc") ;
	
	my %seed_ps; my %seedseq_ps;
	while (<FPS>) {
		my @t = split (/\t/, $_);
		if ($t[0] eq "H") {
			chomp ($t[9]); chomp ($t[8]);
			$seed_ps{$t[9]}++ ;
			$seedseq_ps{$t[9]}.= $t[8].", " ;
		}
	
		if ($t[0] eq "S") {
			$seed_ps{$t[8]}++ ;
			$seedseq_ps{$t[8]}.= $t[8].", ";
		}
	}
	close FPS;
    
	print RPS "OTU\tSeed sequence\t#Sequences\tSequences\n";
	my $ips=1;
	foreach my $e ((sort { $seed_ps {$b} <=> $seed_ps {$a}} keys %seed_ps) ){
		print RPS "$ips\t$e\t$seed_ps{$e}\t$seedseq_ps{$e}\n";
		$ips++;
		print SEQPS ">$e\n$seq{$e}\n";
	}
	
	close RPS;
	close SEQPS;
	
#############################"
	print "\nSamples description printed in index_$Clst file.\nOTUs of environments clustered together are printed in $pooled_sample/pooled_sample_SEQ_OTU file.\nYou can now process $pooled_sample/pooled_sample_SEQ_OTU with panam.pl\n\n";

}

